struct Linbcg {
	virtual void asolve(VecDoub_I &b, VecDoub_O &x, const Int itrnsp) = 0;
	virtual void atimes(VecDoub_I &x, VecDoub_O &r, const Int itrnsp) = 0;
	void solve(VecDoub_I &b, VecDoub_IO &x, const Int itol, const Doub tol,
		const Int itmax, Int &iter, Doub &err)
	{
		Doub ak,akden,bk,bkden=1.0,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
		const Doub EPS=1.0e-14;
		Int j,n=b.size();
		VecDoub p(n),pp(n),r(n),rr(n),z(n),zz(n);
		iter=0;
		atimes(x,r,0);
		for (j=0;j<n;j++) {
			r[j]=b[j]-r[j];
			rr[j]=r[j];
		}
		//atimes(r,rr,0);
		if (itol == 1) {
			bnrm=snrm(b,itol);
			asolve(r,z,0);
		}
		else if (itol == 2) {
			asolve(b,z,0);
			bnrm=snrm(z,itol);
			asolve(r,z,0);
		}
		else if (itol == 3 || itol == 4) {
			asolve(b,z,0);
			bnrm=snrm(z,itol);
			asolve(r,z,0);
			znrm=snrm(z,itol);
		} else throw("illegal itol in linbcg");
		while (iter < itmax) {
			++iter;
			asolve(rr,zz,1);
			for (bknum=0.0,j=0;j<n;j++) bknum += z[j]*rr[j];
			if (iter == 1) {
				for (j=0;j<n;j++) {
					p[j]=z[j];
					pp[j]=zz[j];
				}
			} else {
				bk=bknum/bkden;
				for (j=0;j<n;j++) {
					p[j]=bk*p[j]+z[j];
					pp[j]=bk*pp[j]+zz[j];
				}
			}
			bkden=bknum;
			atimes(p,z,0);
			for (akden=0.0,j=0;j<n;j++) akden += z[j]*pp[j];
			ak=bknum/akden;
			atimes(pp,zz,1);
			for (j=0;j<n;j++) {
				x[j] += ak*p[j];
				r[j] -= ak*z[j];
				rr[j] -= ak*zz[j];
			}
			asolve(r,z,0);
			if (itol == 1)
				err=snrm(r,itol)/bnrm;
			else if (itol == 2)
				err=snrm(z,itol)/bnrm;
			else if (itol == 3 || itol == 4) {
				zm1nrm=znrm;
				znrm=snrm(z,itol);
				if (abs(zm1nrm-znrm) > EPS*znrm) {
					dxnrm=abs(ak)*snrm(p,itol);
					err=znrm/abs(zm1nrm-znrm)*dxnrm;
				} else {
					err=znrm/bnrm;
					continue;
				}
				xnrm=snrm(x,itol);
				if (err <= 0.5*xnrm) err /= xnrm;
				else {
					err=znrm/bnrm;
					continue;
				}
			}
			if (err <= tol) break;
		}
	}

	Doub snrm(VecDoub_I &sx, const Int itol)
	{
		Int i,isamax,n=sx.size();
		Doub ans;
		if (itol <= 3) {
			ans = 0.0;
			for (i=0;i<n;i++) ans += SQR(sx[i]);
			return sqrt(ans);
		} else {
			isamax=0;
			for (i=0;i<n;i++) {
				if (abs(sx[i]) > abs(sx[isamax])) isamax=i;
			}
			return abs(sx[isamax]);
		}
	}
};